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Q-i Stochastic differential equations, especially the one called Langevin equation, 

Qh play an important role in many fields of modern science. In this paper, we 

use the bicolour rooted tree method, which is based on the stochastic Tay- 
lor expansion, to get the systematic pattern of the high order algorithm for 
Langevin equation. We propose a popular test problem, which is related 
O to the energy relaxation in the double well, to test the validity of our al- 

gorithm and compare our algorithm with other usually used algorithms in 
ph simulations. And we also consider the time-dependent Langevin equation 

Pi with the Ornstein-Uhlenbeck noise as our second example to demonstrate 

the versatility of our method. 
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1. Introduction 

O 

Nature is full of randomness, from nucleus to whole galaxy, from inor- 
ganism to organism and from the domain of science and technology to the 
human society [1, 2, 3, 4, 5, 6, 7, 8]. Although the mechanisms of randomness 
^ are different from one field to another, the ways to describe them are similar. 

The stochastic differential equation (SDE) is a good approach to describe 
the randomness. The earliest work on SDEs was done to describe Brownian 
motion in Einstein's famous paper and at the same time by Smoluchowski. 
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Later Ito and Stratonovich put SDEs on a more solid theoretical foundation. 
In 1908, a French physicist, Paul Langevin, developed an equation called the 
Langevin equation (LE) thereafter, which incorporated a random force into 
the Newton equation, to describe the Brownian motion. Langevin equation 
is an equation to mechanics using simplified models and using SDEs to ac- 
count for omitted degrees of freedom. There are many branches with rich 
contents which have been derived in the last 100 years. For example, the 
reaction kinetic dynamics in chemistry [9], the molecular motor and protein 
folding in biology [10, 11], the intracellular and intercellular calcium signal- 
ing, quantum Brownian motion and the stochastic quantization in physics 
[12, 13, 14], even the stock market fluctuations and crashes [8] are all related 
to the Langevin equation. The Langevin equation plays an important role in 
modern science, however only a few of them can be analytically solved, thus 
it is necessary to develop a numerical algorithm which incorporates both the 
computation efficiency and accuracy. 

The general structure of the stochastic differential equation is 

x i = f i (X(t))+g i {X(t))Z i (t) (z = l,2,...,rf), (1) 

where X(t) = {x±, X2, xj}, fi(X(t))s are the deterministic part of the 
equations of motion, gi(X(t))s are the diffusion coefficients and £i(t)s are a 
set of independent gaussian random variables with correlation function 

<&(*)&(*') >=SijS(t-t). (2) 

To get a certain order algorithm for the SDE, we can directly do the 
stochastic Taylor expansion of Eq.(l) to our desired accuracy [15, 16]. This 
method is very explicit and suits for many cases of the SDEs, however, since 
this expansion is too laborious to generate to high orders, we need to find 
a systematic pattern to overcome such difficulty In this paper, we use the 
bicolour rooted tree method (BRT) based on the stochastic Taylor expansion 
to obtain the high order algorithm for SDEs systematically. 

In the field of numerical method for solving ordinary differential equa- 
tions, J. C. Butcher develops a rooted tree method which relates each term 
in the ordinary Taylor expansion to a rooted tree [17]. His method can excel- 
lently make the laborious ordinary Taylor expansion systematic in a heuristic 
way. Then K. Burrage and P. M. Burrage expand the rooted-tree method 
to the bicolour rooted tree method which relates each term in the stochastic 
Taylor expansion to a bicolour rooted tree [16] for the sake of solving SDE. 
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They give an explicit Runge-Kutta method of order 2.5 in their paper for 
the SDE. In this paper, we further develop their works, propose a new type 
of the bicolour rooted tree method, and apply it to the LE case. 

Since the intricacy of the numerical method for SDE, the order of it is 
heretofore not great than 2.5 [16, 18]. But for some special kinds of the SDE, 
for example, the Langevin equation, a high order algorithm can be acquired. 
Hershkovitz has developed a fourth order algorithm for the LE [15], which 
is based on the stochastic Taylor expansion. In this paper, we use the BRT 
method to improve the accuracy to order 7 of deterministic part and order 
4.5 of stochastic part (o (7,4.5)). 

In section 2, we briefly introduce the BRT method and explore the rela- 
tion between the terms in the stochastic Taylor expansion and the bicolour 
rooted trees. We find that the stochastic Taylor expansion is just equal to 
the sum of all the non-isomorphic bicolour rooted trees. In section 3, due 
to the structure of LE, we can use the BRT method to obtain our algo- 
rithm o (7,4.5) for the LE. In section 4, we use two examples to verify the 
validity and demonstrate the versatility of our algorithm. The first one is 
the energy relaxation in the double well. We compare our results with the 
previous results obtained by other algorithms and show the convergence of 
these different algorithms. The second one we present an algorithm for the 
time-dependent Langevin equation with the Ornstein-Uhlenbeck noise, and 
our results are readily agreed with the previous ones. 

2. Bicolour rooted tree method 

To cope with the intricacy of the Taylor expansion of SDE, a method 
which is called bicolour rooted tree method (BRT) [16] based on the rooted 
tree method [17] developed by J. C. Butcher is introduced to conveniently 
do the stochastic Taylor expansion of SDE. 

Let us firstly transform the Eq.(l) into the following equation, 



where Wi(t) is the Wiener process, and the symbol o implies that the SDE 
considered in this paper is in the Stratonovich sense, for the Stratonovich 
integral satisfies the usual rules of calculus [18]. One can therefore integrate 
Eq.(3) from to h, 



dxi = fi(X(t))dt + 9i(X(t)) o dW.it) (i = 1, 2, .., d), 



(3) 




(4) 
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Taylor expansion of the functions gives, 

d 



/« = E htit ^m^-rmm, (5) 



n=0 m=l 
oo d 

*(*(*)) = £ -y(E ^ «— ) B ft(*(0)). (6) 

n=0 m=l 

Now taking the last two equations into Eq.(4), one can easily get 

8xi(h) = Xi(h) - Xi(Q) = 
h 1 



(P + f}Sxj(s) + yfjk Sx j( s ) Sx k(s) + ■■■ )ds+ ^ 
1 . 1 . 

(g l + g)8xj(s) + -g l jk Sx j (s)5x k (s) + • • • ) o dWi(s), 

where p = fi(X(0)), fj k = ^:^-fi(X(0)) and the repeated indices except 
i (the number of the equations) imply the Einstein's summation convention 
throughout the paper. 

Then the terms with 0th derivative in Eq.(7) are, 

6x° i (h) = f i J , h + g i J ith , (8) 

where J ,/i = Jo ds and = J Q h odWi(s), so 5xi(h) = 5x®(h) + • • ■ , substi- 
tuting it for Eq.(7) gives the 1st derivative terms, 

Sxl(h) = /;/ ! •/„„,, + fig j J j0>h + g)f j J oi>h + gtfJfrh, (9) 

where Jj 1 j 2 -j k ,t is the Stratonovich multiple integral [18], and the integration 
is with respect to ds if ji — or odWi(s) if ji = i, for example, 

J 012 t = [ odW 2 { Sl ) f odW x {s 2 ) [ 2 ds 3 . (10) 
Jo Jo Jo 

Replacing Eq.(7) by 5xi(h) = Sx^(h) + Sxj(h) + ■ ■ ■ , one can get the 2nd 
derivative terms Sxf(h) and performing this procedure recursively will gener- 
ate all the derivative terms in principle. However, close calculation of these 
derivative terms reveals that the complexity will increase drastically as the 
order rises. For this reason, we adopt the BRT method developed by J. C. 
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Butcher and P. M. Burrage to express each derivative term systematically 
and graphically. 

We will first introduce some useful notations [16]. Take the bicolour 
rooted tree t in Fig.l as an example, The tree has 8 vertices, each vertex can 
be colored by white node (o) or black node (•) which is the representative 
of stochastic node (a) or deterministic node (r). If ti, ■ ■ • ,t m are bicolour 
rooted trees, then [t\, ■ ■ ■ , t m ] and {t±, • • • , t m } are trees in which ti, ■ ■ ■ ,t m 
are each joined by a single branch to • or o, respectively. We can therefore 
rewrite the tree t in a compact form [a, [{cr}, [r, a}}}. To conveniently calculate 
the weight of this tree, we define the following terms: the degree of the vertex 
C(v)(v G t) in the BRT is equivalent to the degree of vertex D{v) in the graph 
theory except the root 1 with C(l) = 1 + D(i). Sis the symmetry factor of 
the tree t, for example, the trees interchanged the branches joint to vertex 1 
or 3 or 5 are regarded as identical with tree t, therefore the symmetry factor 
of tree t is 2x2x2 = 8. Then tracing the stochastic Taylor expansion of 
the Eq.(7), we find that the elementary weight of the tree, which is also the 
coefficient of each term in the expansion, is 

n 

where n is the total number of vertex in the tree and vertex Vi € t. Now we 
introduce the elementary derivative and elementary integral here [16]. An 
elementary derivative F(t) can be associated with a BRT such that 

F(r) = f,F(a)=g, 

' f( m )[F(ti), ■ ■ ■ ,F(t m )],t— [ti,- ■ ■ ,t m ] (12) 
gW[F(ti),--- ,F{t m )],t = {t u .-. ,t m }, 

and the elementary integral can be written as 

Qs(t) = Jo,s,Q s {cr) = J ijS , 



F(t) 



f°a!u(Y[ 9 u (tj)),t= [ti, ■■■ ,t r 



(13) 



9 3 (t) - 

Fig. 1(a) illustrates the elementary weight, derivative and integral graphically 



fo°dw i (u)(ue u (t j )),t= {ti, ■ ■ ■ ,t m }. 
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8 t = [a,[{a},[r,a}]] 

5 = * *(fi) 8 (fi)^)' = l 

F{t) = f"(9,f"Ql9,f"{f,9))) 

= fik9 j fL9 l n9 n fn P f n 9 P 
8h(t) = J dsi,/ JjSl / ds 2 J r a,s 2 I dszJa, S3 J v>! 











(a) 



8 *,.(//) = • + o + 
+ ■ ■ ■ 



+ 



+ 



• o 

+ 



4 4 6 6 



(b) 

Figure 1: Illustration of the BRT method. Fig. 1(a) shows a bicolour rooted tree t, and its 
elementary weight a(t), derivative F(t) and integral 9 h (t), respectively. Fig. 1(b) shows the 
Oth and 1st derivative terms of stochastic Taylor expansion of Sxi(h) by the BRT method. 
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Therefore the stochastic Taylor expansion is given by 
6xi(h)=Y,a(t)Ht)0h(t) 

=sum of all the non-isomorphic 
bicolour rooted trees, 

where T is the set of non-isomorphic bicolour rooted trees. Fig. 1(b) illustrates 
how to use this formula to express Sxi(h) = Sx^(h) + Sx](h) + • ■ • . 

3. Algorithm 

Due to the complexity of stochastic Taylor expansion, we only consider 
Langevin equation (LE) which plays an important part in the fields involving 
randomness in this paper. An N dimensional set of coupled LEs has the form 

,_ -^ + 6(0, (18) 



11 o 

oqi 

where V(q(i)) is the external potential, %(i — 1, • • • , N) is a set of friction 
parameters, is random noise with zero mean, the correlation relation is 

<&(*)&(*') >=2 li TS ij S{t-t'), (16) 
and the Hamilton canonical equations are 

(it) 

The form of Eq.(17) where only every second equation has a noise term with 
constant diffusion coefficient, as well as the potential V(q(i)) is unrelated to 
p(t), makes it possible to sharply decrease the complexity of Eq.(14) so as 
to obtain a high order algorithm for LE. 

For the Eq.(17), we can translate it into the form 

x i = MX(t))+g i Z i (t) (i = 1,2, ...,27V) 
<ti(t)ti(t')>=5 ij 5(t-t'), 

where fi(X(t)) is equal to x i+ i for odd % and to — dV(X.(t))/dxi_i — jiXi 
for even i, X(£) = {x±, X3, X2N-1}, 9i is a set of constants with ^ = 
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if i is odd number. Because of the property of V(q{t)), one can find that 
fjk 7^ on ly ^ j an d & are both odd numbers. From above properties, a 
key property for the simplification of the stochastic Taylor expansion, that 
is, ■ • • fjkdj ■ ■ • = 0, can be found. We can rewrite it in the compact form as 
follow: 

•••,[••-, a,- ••],-•• = 0, (19) 

so if a bicolour rooted tree has this structure, it should have no contribution 
to the stochastic Taylor expansion. 

From the analysis above, one can obtain a general method for solving the 
Langevin equation numerically. If we want to get a numerical method to the 
order o(m,n), we should: 

(a) For the deterministic part: 

Solve it by the standard Runge-Kutta method of order m. 

(b) For the stochastic part: 

(i) Write down all the non-isomorphic bicolour rooted trees that can avoid 
the structure (19); 

(ii) Attach each vertex with white or black so as to make the tree have 
order n. 

(c) Add up the results of (a) and (b). 

Using these three criteria, all the terms up to order o (7,4.5) are, 

6 

5 Xidet (h)=^2RK{j), 
Sx iran (h) =a + [a] + [[a]} + [[[a]}} + [r, [a]} 

+ [[[[*]]]] + [t,[H] + [MM 

+ [r,r, M] + [[t,[<7]]] + H[4 
Sxi(h) =5x idet {h) + 8x iran (h). 

where the RK(j), (j = 0, • • • ,6) are the BRTs with j + 1 deterministic nodes 
only which are identical to the terms of standard Runge-Kutta method for 
ODEs. These terms compose the deterministic part of our algorithm, and 
we can use Runge-Kutta method to solve the deterministic part numerically 
[17,19]. 

Then we try to find a way to calculate the stochastic part in Eq.(20). 
We here introduce a method to approximate the elementary integral which 
is developed by P. E. Kloeden and E. Platen [18]. They showed that if 
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(20) 



W(£) = (Wi(t), • • • , W N (t)) is an N-dimensional Wiener process on the time 
interval [0, h], the componentwise Fourier expansion of the Brownian bridge 
process W(t) - |W(/i) is 



wt(t) - jWtih) = -a l0 + cos ("tt) + % sin (^r))' ( 21 ) 



where 0^,6^ are N(0, h/27r 2 j 2 ) distributed and pairwise independent, then 



setting t = in the equation (21) gives = — 2 = a° 



Now, we begin to calculate the stochastic part of Eq.(20). Firstly, let us 
set the W(0) = 0, then use equation (21), we can easily find that 



a = a(a)F(a)6 h (a) = g% h = g l W t {h) = g'u], 



M = fi9 j Jjo,h = f}9 3 \ d Sl / odWfa) 



where Wi(h) = Wi is a set of independent Gaussian random variables sampled 
from iV(0, h). Similarly calculation gives 





= irJiriu'4 






[r, H 


= h 3 f jk fft9 l oof 


E?]]\] 


= h 4 "fjfkfl 'fm9 m ^m 


[r, [H] 


= h%f j f?f l m g m u; 7 m 


[[r],H 


= h 4 fjkflf l fm9 muj m 




h 4 

ri fjfkfl m 9 
— 2 J i klJ J Jrn ^ " 








h 3 

— 2 hkh 9 J my lL lm 




(22) 



hf]g J {- 
W<4 



Wj(h) 
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where 



4 _ Wi a° a* 6J 



6 4 



* 120 48 8 12 8 
7 a° a] bj ft? 

w ^30 + ii + ¥ + ¥ + ¥ 



W^ a^_a[ 6[_6f (24) 
20 16 8 12 8 



a; 9 = — - + — + — 

10 8 4 6 4 



Wi a° 6J 6? 
—+—+—+— 
40 24 12 4 

a° = -2^a^-, &} = bij, 

3=1 3=1 



E % = E (-j^ 6 - 



and the first non- Gaussian random variable is 

i ,,i,,9, n,5 W M a °S° a M, 6 i 6 i 
^■-^■+0,0;,.-— ^ + ^ 



oo 



— E 1 i,\2^ aik<1 3 k + bikbjk)- 

k=l 



(25) 



We can find that there are only 5 independent variables among uj, (j = 



10 



1, • • • ,10). Let's choose oj],u 2 , uif, u>f, uf as the independent variables, then 



wj = - 

"f = + (26) 

^ = w? - 2uf + 2wf 
^° = ^ - 2cf 

Now, the last procedure we should do is to determine the five Gaussian 
random variables Wi,a®, a], b], b 2 and the non-Gaussian random variable c^-. 
We truncate Cij to the first term, that is, rs (aa%i + hibji) /n 2 ■ Since 
Wi,dij,bij are independent and Wj ~ N(0,h), aij ~ N(0, h/2n 2 j 2 ), bij ~ 
iV(0, h/2n 2 j 2 ), we can see that 



(27) 



a, 1 ~M0, — ), 6, 2 ~ MO, — ^— ), 
v ' 1890 ; ' * V ' 18900 ; ' 

ni ^ n ^ i ^ 

< ai<»i >= -qq, < a^aa >= < c^au >= — , 

< W >= i^o> < ^ >= ^3, < &&i >= ^5. 

Let 4>\i 4>1i 4>h 4>ii 4>i ) 0j to be seven independent standard Gaussian 
random variables, then use Eq.(27), we can get 

Wi = Vh<t>] 



a 



3 rt 

1 ~ 30 1 v/3 + 



b) 



h 4 (28) 
180^ 



2= 0^ 0f 
°* 63 l v^ + 10 J 

a rt = ^(-0.1754930^ + 0.1393480^ + 0.02109060°) 
b a = ^(0.216350^ + O.O6179950, 5 + 0.005843420?) 
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The BRT method gives an algorithm for the Langevin equation so long as 
we determine the deterministic and the stochastic part of Eq. (20) respectively 
and add up each other. The deterministic part can be solved by the standard 
Runge-Kutta algorithm [19], and the stochastic part can be solved by the 
Eqs.(22)-(23). 

4. Numerical simulations 

4-1. Energy relaxation in double well 

To verify the validity of our algorithm, the Kramers equation will be 
considered as the severe test for our algorithm. The form is as follow: 

p(t) = -v'(i(t))-v(t)+m, 

and is the Gaussian random force obeying the fluctuation dissipation 
theorem 

< mat) >= 2 7 T. (30) 
Our method implies that the algorithm for Eqs.(29)-(30) is, 

q(t + h) =q de t(q(t),p(t), h) + q ran (q(t) , p(t) , h), ^ 
p(t + h) =Pdet(q(t),p(t),h)+p ran (q(t),p(t),h), 

where q<iet(q(t) , p(t) , h) and Pdet(q(t),p(t),h) are the results of evolving the 
equations in the period 0-h by the seventh-order Runge-Kutta algorithm [19] 
which is used in the ODE, and the stochastic part of Eq.(31) is, 

q ran {q{t),p{t), h) =y/¥fr[{h4 - h 2 ju 3 + /i 3 7 V[ - h^ul) 
+ (-h 3 uj* + 2h i 1 ut)V" - h A p(t)uj™V'"l 

Pran(q(t), P (t), h) = ^f[(iU l 2 - hjU^ + h^UJ 3 - h 3 ^ 

+h ^co!) + (-h 2 cu 3 + 2h 3 1 u 4 2 - Zh^cu^V" (32) 
+ (-h 3 p(t)u 2 ) + h 4 -fp(t)u 7 2 + h^p{t)ujl 
+h 4 1 p{t)u\ Q )V'" + h 4 V"v"u 6 2 + h 4 V'v'"co s 2 

-^v"" P 2 (t)u!\-h 3 jTv'"nl 2 , 

where u 2 (i = !,••• , 10) and Vl\ 2 have been defined in the previous section. 
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The double well potential in this example is, 



V{q)=q*-2q\ 



(33) 



It has two minima located at q = ±1 and a potential barrier with the height 
AV^ = 1 between the two wells. The friction coefficient 7 is set to 1. The 
initial condition is chosen on the top of the barrier. The average is taken 
over 5000 realizations of the Gaussian random force during the trajectory. 
Fig. 2 shows the result which is compared with the Euler method and the 
Heun method [20]. We perform these three methods at T=0.05 and T=0.2 
respectively. We find that the results of these three different methods are 
almost agreed. Nevertheless, the step size of our method, Heun method and 
Euler method are 0.1, 0.001 and 0.0001, respectively. The Kramers equation 
has been simulated extensively by many authors (Ref. [15] and the references 
therein). As for the convergence, we compare our algorithm with previous 
algorithms here. Fig. (3) shows the convergence of the three algorithms for 
solving the Kramers equation. It is evident that our algorithm diverges slowly 
than the other algorithms as the step size increases. 

4-2. Stochastic resonance 

Stochastic resonance, which is originally developed to explain the ice ages 
[21, 22], has spread well beyond physics and left its fingerprints in many other 
research areas [23, 24], such as complex networks [25], biological systems [26], 
neuroscience [27, 28] and quantum computing [29] . The governing equations 
in these very different fields are essentially Langevin equation or its general- 
izations. We present an example of stochastic resonance in neuroscience to 
demonstrate our algorithm in the case of time-dependent Langevin equation 
with the Ornstein-Uhlenbeck noise [23]. An enlightening model in the neu- 
ronal dynamical systems is the noise-driven bistable system whose equations 
can be described as follow: 



where e(t) is the Ornstein-Uhlenbeck noise with intensity D and the inverse 
of characteristic time A, and the system is driven by an external periodic 
force with amplitude A and frequency u. 



x + 7± = — V (x) + A cos(ut) + e(t) 
<e(t)e(t) >=DAexp(-A|f-t'|), 



(34) 
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Figure 2: Comparison of different numerical methods for solving the energy relaxation in 
the double well potential. The friction coefficient 7 is equal to 1. The average is taken 
over 5000 realizations for all the algorithms. The solid line is our method with step size 
0.1, the dash line is the Heun method with step size 0.001, and the short dash line is the 
Euler method with the step size 0.0001. Panel a shows the relaxation at high temperature 
T=0.2 and panel b shows the relaxation at low temperature T=0.05, respectively. 
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Figure 3: The convergence of different numerical methods for solving the Kramers equa- 
tion. The friction coefficient 7 is equal to 1. E is the average energy at a certain temper- 
ature and h is the step size. Panel a shows the convergence at high temperature T=0.2 
and panel b shows the convergence at low temperature T=0.05, respectively. 
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To use our algorithm to solve Eq.(34) numerically, we should first trans- 
form it into, 

y = i, 

x = p, 

P = ~1V - V' ( x ) + A cos(uy) + e, (35) 
e = -Xe + \V2D^(t), 

<mat) >=5(t-t). 

Let y — > x±, x — >■ x 2 ,p — > x 3 ,e — > x 4 , we can further simplify Eq.(35) into a 
compact form, 



x i = f i (X(t))+g£ i (t) (i = 1,2,3,4), 

<amt)>=s(t-t), 



(36) 



with 



/i(X(t)) = 1, / 2 (X(t)) = x 3 , U(X(t)) = -Ax 4 , 

/ 3 (X(t)) = -7x3 -y'(ar 2 ) + Acos(wxi) + x 4 , (37) 

9i = 92= 93 = 0, 5-4 = Av 7 ^, 

then one can easily find that property (19) is held again. 

Accordingly, the numerical method of Eqs.(36)-(37) can be written as 
follow: 

xi(t + h) =x ldet (X(t), h) + x lran (X(t), h), 

X 2 (t + K) =X 2det (X(t), h) + X 2ran (X(t), h), 

x?,{t + h) =x 3det (X(t), h) + x 3ran (X(t), h), 
x A (t + h) =x 4 det(X(t), h) + x 4ran (X(t), h), 

where the deterministic part of Eq. (38) accords with the Ronge-Kutta algo- 
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rithm for the ODEs, and the stochastic part of Eq.(38) is 



x lran (X(t),h) = 0, 

X 2r an{X{t), h) = \V2D[{h 2 Uj\ - /* 3 ( 7 + 

+ /i 4 ( 7 2 + 7 A + A 2 H 6 ) - rfulV"], 

x 3ran (X(t), h) = \V2D[(hcol - h 2 ^ + AH 3 

+ /i 3 ( 7 2 + 7 A + A 2 K 4 - /i 4 ( 7 3 + 7 2 A (39) 

+ 7A 2 + A 3 K) + (-h 3 cot + /i 4 (2 7 + X)cof)V" 
-h 4 x 3 (t)ulv'"}, 
x 4ran (X(t), h) = \\f2D[u\ - hXul + h 2 \ 2 u\ 
-h 3 X 3 ut + h 4 X 4 ut}. 

The double well potential in this example is, 

V(x) = - l -x 2 + - A x\ (40) 

and the periodic driven force's amplitude A and frequency u are 0.03 and 
0.01 respectively. We consider the relation between the amplitude of output 
of the system < x > and the noise intensity D. The average is taken over 
5 x 10 5 realizations and the Heun method is used as a comparison. We then 
compare our results with the model mentioned in [23]: 

x = -V'(x)+AcoB(ut)+Z(t), 
<t(t)t(t')>=2D5{t-t'). 

The theoretical result of < x > in this model is < x >= 4 ; 2 n fc =1 where 

r fc = ^exp(-^f). 

Fig. (4) shows the results of our simulations. The black line and the red 
line are the simulations of our method with step size 0.1 and the Heun method 
with step size 0.01 respectively, with the parameters 7 and A equal to 1. We 
now see that the results of our method and the Heun method are almost the 
same, however, the step size of our method is larger than the one used in 
Heun method. The parameters of the green line is the same as the black 
line except A = 10, that is, the characteristic time is shorter, and in this 
condition, the Ornstein-Uhlenbeck noise is closer to the Gaussian noise. We 
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Figure 4: Comparison of different numerical methods for the stochastic resonance in the 
noise-driven bistable system. The periodic driven force's amplitude A and frequency uj 
are 0.03 and 0.01, respectively. The friction coefficient 7 is equal to 1. The amplitude of 
output is averaged over 5 x 10 5 trajectories. The black line is our method with step size 
0.1 and the red line is the Heun method with step size 0.01. Both the two lines have the 
characteristic time 1. The green line is our method with step size 0.1 and characteristic 
time 0.1. The blue line is the theoretical result of Eq.(41). 
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can find that the resonant noisy intensity (the maximum of the line) shifts 
left when we shorten the characteristic time. In other words, lengthening 
the characteristic time can enhance the noise resistance of the system. The 
blue line is the theoretical result of Eq.(41). Since the influence of the inertia 
term x, we can see that the amplitude of output < x > of the theoretical 
result is greater than our numerical result as shown in the green line. 

5. Conclusion 

We have proposed the bicolour rooted tree method to do the stochas- 
tic Taylor expansion systematically. This method can be used to solve the 
stochastic differential equation numerically. In this paper, we focus on the 
Langevin equation which is widely used in modern science. A high order 
algorithm o (7,4.5) is derived in this paper. Comparing with other usual 
algorithms, our method is advantageous in computational efficiency and ac- 
curacy. We present our method in the two examples. In the first example 
of energy relaxation in the double well, we show our method gives the same 
results as presented in other papers, and the convergence is better than the 
other algorithms. In the second example, we propose an algorithm for the 
time-dependent Langevin equation with the Ornstein-Uhlenbeck noise, and 
the result of our algorithm is the same as the one obtained by Heun method. 
It shows our algorithm is suitable for the Langevin equation regardless of the 
time-dependence of the equation. However, the readers should note that we 
only provide the algorithm for Eq.(17) which satisfies the property (19) since 
this property can drastically reduce the complexity of Eq.(14). For the other 
type of SDEs that the property (19) can not be held, such as the Hodgkin- 
Huxley model in neuroscience [27], interested readers can design their own 
algorithms based on the Eq.(14). For the purpose of using our method in 
the more difficult situations, one can consider the case that the diffusion co- 
efficients are variable. All in all, we have provided a systematic scheme for 
searching the high order algorithm for the SDE and find that it can reduce 
drastically when deal with the Langevin equation. 
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